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Abstract 

This paper presents two sufficient conditions to ensure a faithful evaluation of 
polynomial in IEEE-754 floating point arithmetic. Faithfulness means that the 
computed value is one of the two floating point neighbours of the exact result; it 
can be satisfied using a more accurate algorithm than the classic Horner scheme. 
One condition here provided is an a priori bound of the polynomial condition 
number derived from the error analysis of the compensated Horner algorithm. The 
second condition is both dynamic and validated to check at the running time the 
faithfulness of a given evaluation. Numerical experiments illustrate the behavior of 
these two conditions and that associated running time over-cost is really interesting. 

Keywords: Polynomial evaluation, faithful rounding, Horner algorithm, compen- 
sated Horner algorithm, floating point arithmetic, IEEE-754 standard. 

1 Introduction 

1.1 Motivation 

Horner's rule is the classic algorithm when evaluating a polynomial p(x). When performed 
in floating point arithmetic this algorithm may suffer from (catastrophic) cancellations 
and so yields a computed value with less exact digits than expected. The relative accuracy 
of the computed value p(x) verifies the well known following inequality, 

\p(x) — p(x)\ , , 

WK . x , v 71 < a(n) cond(p, x) u. 1 
\p{x)\ 

In the right-hand side of this accuracy bound, u is the computing precision and 
a(n) ~ 2n for a polynomial of degree n. The condition number cond(p, x) that only 
depends on x and on p coefficients will be explicited further. The product a(n) cond(p, x) 
may be arbitrarily larger than 1/u when cancellations appear, i.e., when evaluating the 
polynomial p at the x entry is ill-conditioned. 
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When the computing precision u is not sufficient to guarantee a desired accuracy, 
several solutions simulating a computation with more bits exist. Priest-like "double- 
double" algorithms are well-known and well-used solutions to simulate twice the IEEE- 
754 double precision [HlEj. The compensated Horner algorithm is a fast alternative to 
"double-double" introduced in [2J — fast means that the compensated algorithm should 
run at least twice as fast as the "double-double" counterpart with the same output 
accuracy. In both cases this accuracy is improved and now verifies 

\p(x) - p(x)\ ^ „„ ,/„ \ 2 



\p(x)\ 



< u + j3(n) cond(jo, x) u , (2) 



with (3(n) ~ 4n 2 . This relation means that the computed value is as accurate as the 
result of the Horner algorithm performed in twice the working precision and then 
rounded to this working precision. 

This bound also tells us that such algorithms may yield a full precision accuracy for 
not too ill-conditioned polynomials, e.g., when (3{n) cond(p, x)u < 1. 

This remark motivates this paper where we consider faithful polynomial evaluation. 
By faithful (rounding) we mean that the computed result p{x) is one of the two floating 
point neighbours of the exact result p(x). Faithful rounding is known to be an interesting 
property since for example it guarantees the correct sign determination of arithmetic 
expressions, e.g., for geometric predicates. 

We first provide an a priori sufficient criterion on the condition number of the poly- 
nomial evaluation to ensure that the compensated Horner algorithm provides a faithful 
rounding of the exact evaluation (Theorem in Section EJ). We also propose a validated 
and dynamic bound to prove at the running time that the computed evaluation is ac- 
tually faithful (Theorem El in Section HJ. We present numerical experiments to show 
that the dynamic bound is sharper than the a priori condition and we measure that the 
corresponding over-cost is reasonable (Section EJ). 



1.2 Notations 

Throughout the paper, we assume a floating point arithmetic adhering to the IEEE-754 
floating point standard [5]. We constraint all the computations to be performed in one 
working precision, with the "round to the nearest" rounding mode. We also assume that 
no overflow nor underflow occurs during the computations. Next notations are standard 
(see [4, chap. 2] for example). F is the set of all normalized floating point numbers and u 
denotes the unit roundoff, that is half the spacing between 1 and the next representable 
floating point value. For IEEE-754 double precision with rounding to the nearest, we 
have u = 2~ 53 rs 1.11 • 10~ 16 . We define the floating point predecessor and successor of a 
real number r as follows, 

pred(r) = max{/ 6 F// < r} and succ(r) = min{/ G F/r < /}. 

A floating point number / is defined to be a faithful rounding of a real number r if 

pred(/) < r < succ(/). 
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The symbols 0, 0, (g) and represent respectively the floating point addition, subtrac- 
tion, multiplication and division. For more complex arithmetic expressions, fl(-) denotes 
the result of a floating point computation where every operation inside the parenthesis is 
performed in the working precision. So we have for example, a © b = fl(a + b). 

When no underflow nor overflow occurs, the following standard model describes the 
accuracy of every considered floating point computation. For two floating point numbers 
a and b and for o in {+, — , x, /}, the floating point evaluation fl(a 06) of a o b is such 
that 

fl(ao&) = (ao6)(H-e x ) = (ao6)/(l + e 2 ),with |ei|,|e 2 |<u. (3) 

To keep track of the (1 + e) factors in next error analysis, we use the classic (1 + 9k) 
and 7^ notations jH chap. 3]. For any positive integer k, 9k denotes a quantity bounded 
according to 

ku 

\Vk\ < 7fc = :— . 

1 — ku 

When using these notations, we always implicitly assume ku < 1. In further error 
analysis, we essentially use the following relations, 

(l + fc )(l + 0j) < (1 + Ok+j), fcu< 7fc , 7 fc <7 fc+ i. 

Next bounds are computable floating point values that will be useful to derive dynamic 
validation in Section HJ We denotes &(jk) = (ku) 0(10 ku) by j k . We know that 
fl(&u) = ku G F, and ku < 1 implies fl(l — ku) = 1 — ku G F. So 7 fc only suffers from a 
rounding error in the division and 

7fc<(l + u)7 fc - (4) 
The next bound comes from the direct application of Relation (jSJ). For x G F and n G N, 

(i +u rw< fl (3— J^L_). ( 5) 

2 From Horner to compensated Horner algorithm 

The compensated Horner algorithm improves the classic Horner iteration computing a 
correcting term to compensate the rounding errors the classic Horner iteration gener- 
ates in floating point arithmetic. Main results about compensated Horner algorithm are 
summarized in this section; see [2] for a complete description. 

2.1 Polynomial evaluation and Horner algorithm 

The classic condition number of the evaluation of p(x) = J2i=o aiX * a ^ a gi ven data x is 

condfn x) = \ ai ^ 1 = M^L (Q) 

For any floating point value x we denote by Horner (p, x) the result of the floating point 
evaluation of the polynomial p at x using next classic Horner algorithm. 
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Algorithm 1. Horner algorithm 
function r = Horner (p, x) 

for i = n — 1 : — 1 : 

Ti = r i+ i ® x®<Xi 
end 

The accuracy of the result of Algorithm verifies introductory inequality (pj with 
a n u = 72n and previous condition number (0). Clearly, the condition number cond(p, x) 
can be arbitrarily large. In particular, when cond(p, x) > l/j2n, we cannot guarantee 
that the computed result Horner (p, x) contains any correct digit. 

We further prove that the error generated by the Horner algorithm is exactly the sum 
of two polynomials with floating point coefficients. The next lemma gives bounds of the 
generated error when evaluating this sum of polynomials applying the Horner algorithm. 

Lemma 1. Let p and q be two polynomials with floating point coefficients, such that 
p( x ) = Sr=o aiX% an d l( x ) = Xir=o biX 1 . We consider the floating point evaluation of (p + 
q)(x) computed with Horner (p © q, x). Then, in case no underflow occurs, the computed 
result satisfies the following forward error bound, 

\(p + q) (x) - Horner (p®q,x)\ < 7 2n +i ( P + q) (x) ■ (7) 

Moreover, if we assume that x and the coefficients of p and q are non-negative floating 
point numbers then 

(p + q)(x) < (1 + u) 2n+1 Horner (p ®q,x). (8) 

Proof. The proof of the error bound (jJJ) is easily adapted from the one of the Horner 
algorithm (see 4, p. 95] for example). To prove (jHJ) we consider Algorithm ^ where 

T n — «n © b n and r, = r i+ i © x © (a^ © bj) for i — n — 1, . . . , 0. 

Next, using the standard model (jHJ) it is easily proved by induction that, for i — 0, . . . , n, 

i 

J^ian-i+j + b n „ i+] )x 3 < (1 + u) 2m r n _i, (9) 
j=o 

which in turn proves (jHJ) for i = n. □ 
2.2 EFT for the elementary operations 

Now we review well known results concerning error free transformation (EFT) of the 
elementary floating point operations +, — and x. 

Let o be an operator in {+, — , x}, a and b be two floating point numbers, and x = 
fl(a o b). Then their exist a floating point value y such that 

a o b = x + y. (10) 
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The difference y between the exact result and the computed result is the rounding error 
generated by the computation of x. Let us emphasize that relation (fTUJ) between four 
floating point values relies on real operators and exact equality, i.e., not on approximate 
floating point counterparts. Ogita et al. [8] name such a transformation an error free 
transformation (EFT). The practical interest of the EFT comes from next Algorithms^ 
and 0] that compute the exact error term y for o = + and o = x . 

For the EFT of the addition we use Algorithm the well known TwoSum algorithm 
by Knuth jB] that requires 6 flop (floating point operations). For the EFT of the product, 
we first need to split the input arguments into two parts. It is done using Algorithm El of 
Dekker £Q where r = 27 for IEEE- 754 double precision. Next, Algorithm 0] by Veltkamp 
(see P]) can be used for the EFT of the product. This algorithm is commonly called 
TwoProd and requires 17 flop. 



Algorithm 2. EFT of the sum of two floating point numbers. 

function [x, y] = TwoSum (a, b) 

x = a Q)b 
z = x a 

y=(aQ(xQz))@ (b z) 
Algorithm 3. Splitting of a floating point number into two parts. 

function [x,y] = Split (a) 
z = a (2 r + 1) 

x = z (z a) 
y = a x 

Algorithm 4. EFT of the product of two floating point numbers. 

function [x, y] = TwoProd (a, b) 
x = a b 
[a h ,ai] = Split (a) 
[6 ft ,6,]=Split(6) 

y = ai®biQ (((x a h b h ) a/ b h ) Qa h ® b{) 



The next theorem exhibits the previously announced properties of TwoSum and 
TwoProd. 

Theorem 2 ([8]). Let a,b in¥ and x,y G F such that [x, y] = TwoSum(a, b) (Algorithm 
H}). Then, ever in the presence of underflow, 

a + b = x + y, x = aQ)b, \y\<u\x\, \y\<u\a + b\. 

Let a, b G F and x, y G F such that [x,y] = TwoProd(a,6) (Algorithm^. Then, if no 
underflow occurs, 

a x b = x + y, x = a b, \y\ < u\x\, \y\ < u\a x b\. 
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We notice that algorithms TwoSum and TwoProd only require well optimizable floating 
point operations. They do not use branches, nor access to the mantissa that can be time- 
consuming. We just mention that significant improvements of these algorithms are defined 
when a Fused-Multiply-and-Add operator is available [2]. 

2.3 An EFT for the Horner algorithm 

As previously mentioned, next EFT for the polynomial evaluation with the Horner algo- 
rithm exhibits the exact rounding error generated by the Horner algorithm together with 
an algorithm to compute it. 

Algorithm 5. EFT for the Horner algorithm 
function [so,Pn,Po-] = EFTHornerQo, x) 

for i — n — 1 : — 1:0 

[pi,7Ti] = TwoProd(s i+ i,x) 

[si,(Ji] = TwoSum(pi,aj) 

Let 7Tj be the coefficient of degree i in p n 

Let <7j be the coefficient of degree i in p a 
end 

Theorem 3 ([2]). Let p(x) = Y^=o a i xl ^ e a polynomial of degree n with floating point 
coefficients, and let x be a floating point value. Then Algorithm^ computes both 

i) the floating point evaluation Horner (p,x) and 

ii) two polynomials p n and p a of degree n — 1 with floating point coefficients, 
such that 

[Horner (p, x) , p n , p a ] = EFTHorner (p, x) . 

If no underflow occurs, 

p(x) = Horner (p,x) + (p* +p a )(x). (11) 

Moreover, 

(Pn +Pa){x) < l2nP{x)- (12) 

Relation (|TT| means that algorithm EFTHorner is an EFT for polynomial evaluation 
with the Horner algorithm. 

Proof of Theorem^ Since TwoProd and TwoSum are EFT from Theorem |21 it follows 

that s i+ ix = Pi + 7Tj and Pi + cij = + <7j. Thus we have Sj = s i+ \X + aj — 7Tj — a iy for 
i — 0, . . . , n — 1. Since s n = a n , at the end of the loop we have 

n n— 1 n— 1 

so = ^2 aiX% ~ 7X1x1 ~ ° iX% ' 

i=0 i=0 i=0 

which proves (fTTj) . 
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Now we prove relation (jl2j) According to the error analysis of the Horner algorithm 
(see |H p. 95]), we can write 

n-1 

Horner (p, x) = (1 + # 2 „)a„£ n + ^(1 + 9 2l+ i)aix\ 

i=0 

where every 9k satisfies \9k\ < 7fc- Then using (fTTj) we have 

n-1 

(Pn + Pa)(x) = p{x) - Horner (p, x) = -^„a„a; n - 2J 02i+ia»z*. 

i=0 

Therefore it yields next expected inequalities between the absolute values, 

n-1 

(Ptt +Pa){x) < l2n\a n \\x\ n + 72i+l | Oi 1 1 Z | * < 72nP(^). 

i=0 

□ 

2.4 Compensated Horner algorithm 

From TheoremOlthe final forward error of the floating point evaluation of p at x according 
to the Horner algorithm is 

c = p(x) - Horner (p,x) = (p w +p a )(x), 

where the two polynomials p w and p a are exactly identified by EFTHorner (Algorithm EJ) 
— this latter also computes Horner (p, x). Therefore, the key of the compensated algorithm 
is to compute, in the working precision, first an approximate c of the final error c and 
then a corrected result 

f = Horner (p, x) © c. 

These two computations leads to next compensated Horner algorithm CompHorner (Al- 
gorithm EJ) . 

Algorithm 6. Compensated Horner algorithm 

function r = CompHorner (p, x) 
[f,p w ,p CT ] = EFTHorner (p, x) 
c = Horner (p n © p a , x) 

f = r © c 

We say that c~is a correcting term for Horner (p, x). The corrected result f is expected 
to be more accurate than the first result Horner (p, x) as proved in next section. 

3 An a priori condition for faithful rounding 

We start proving the accuracy behavior of the compensated Horner algorithm we pre- 
viously mentioned with introductory inequality (j2J) and that motivates the search for a 
faithful polynomial evaluation. This bound (and its proof) is the first step towards the 
proposed a priori sufficient condition for a faithful rounding with compensated Horner 
algorithm. 
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3.1 Accuracy of the compensated Horner algorithm 

Next result proves that the result of a polynomial evaluation computed with the compen- 
sated Horner algorithm (Algorithm |UJ) is as accurate as if computed by the classic Horner 
algorithm using twice the working precision and then rounded to the working precision. 

Theorem 4 (|2J). Consider a polynomial p of degree n with floating point coefficients, 
and x a floating point value. If no underflow occurs, 

|CompHorner (p, x) — p(x)\ < u\p(x)\ +72„p(^)- (13) 

Proof. The absolute forward error generated by Algorithm El is 

| r — p(x)\ — |(r © c) — p(x)\ — |(1 + e)(r + c) — p(x)\ with |e| < u. 

Let c = (p n + p a )(x). From Theorem |3] we have r = Horner (p, x) = p(x) — c, thus 

| r — p(x)\ = |(1 + e) (p(x) — c + c) — p(x)\ < u\p(x)\ + (1 + u)| c — c|. 

Since c= Horner (p n (B p a , x) with p n and p a two polynomials of degree n — 1, Lemma ^ 
yields |c- c| < ^n-iiPn + Pa)(x). Then using (JUJ) we have \c- c\ < 7 2n -i72nP r ( a; )- 
Since (1 + u)7 2n _i < 72 n , we finally write the expected error bound (|T3|) . □ 

Remark 1. For later use, we notice that | c — c\ < 72n-i72nP(x) implies 

I c — c| <-y% n P(x). (14) 

It is interesting to interpret the previous theorem in terms of the condition number 
of the polynomial evaluation of p at x. Combining the error bound ()13|) with the con- 
dition number © of polynomial evaluation gives the precise writing of our introductory 
inequality (J2J), 

|CompHorner Qo, x) — p(x)\ ^ 2 



x , rj < u + 7 2n cond(p,x). (15) 

In other words, the bound for the relative error of the computed result is essentially 
times the condition number of the polynomial evaluation, plus the inevitable summand 
u for rounding the result to the working precision. In particular, if cond(p, x) < 
then the relative accuracy of the result is bounded by a constant of the order u. This 
means that the compensated Horner algorithm computes an evaluation accurate to the 
last few bits as long as the condition number is smaller than u/7^ w l/4n 2 u. Besides 
that, relation (fTKl) tells us that the computed result is as accurate as if computed by 
the classic Horner algorithm with twice the working precision, and then rounded to the 
working precision. 

3.2 An a priori condition for faithful rounding 

Now we propose a sufficient condition on cond(p, x) to ensure that the corrected result 
r computed with the compensated Horner algorithm is a faithful rounding of the exact 
result p(x). For this purpose, we use the following lemma from [TP] . 
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Lemma 5 ( [10J ) . Let r,5 be two real numbers and r = fl(r). We assume here that r is 
a normalized floating point number. If \5\ < ^ \ r| then r is a faithful rounding ofr + S. 

From Lemma El we derive a useful criterion to ensure that the compensated result 
provided by CompHorner is faithfully rounded to the working precision. 

Lemma 6. Let p be a polynomial of degree n with floating point coefficients, and 
x be a floating point value. We consider the approximate r of p(x) computed with 
CompHorner (p, x), and we assume that no underflow occurs during the computation. Let 
c denotes c = (p n +p a )(x). If\c — c\ < f |r|, then r is a faithful rounding ofp(x). 

Proof. We assume that \c — c| < j\f\. From the notations of Algorithm El we recall 
that fl( r + c) = r. Then from Lemma El it follows that r is a faithful rounding of 
r+c + c—c= r + c. Since [r,p n ,p a ] = EFTHorner (p, x), Theorem 01 yields p(x) = r + c. 
Therefore r is a faithful rounding of p(x). □ 

The criterion proposed in Lemma El concerns the accuracy of the correcting term c. 
Nevertheless Relation pointed after the proof of Theorem H] says that the absolute 
error \c — c\ is bounded by ^i\ n p[x). This provides us a more useful criterion, since 
it relies on the condition number condQo, x), to ensure that CompHorner computes a 
faithfully rounded result. 

Theorem 7. Let p be a polynomial of degree n with floating point coefficients, and x a 
floating point value. If 

cond(jP, x) < - - " u7 2n ~ 2 , (16) 

then CompHorner (p, x) computes a faithful rounding of the exact p(x). 

Proof. We assume that (|T6*|) is satisfied and we use the same notations as in Lemma El 

First we notice that r and p(x) are of the same sign. Indeed, from ()13|) it follows 
that \ f/p(x) — 1| < u + 7| n cond(p, x), and therefore r/p(x) > 1 — u — cond(p, x). 
But (fTT)|) implies that 1 — u — 7| n cond(p, x) > 1 — 3u/(2 + u) > 0, hence f/p(x) > 0. 
Since r and p(x) have the same sign, it is easy to see that 

(1 -u)\p(x)\ -ll n p{x) < \f\. (17) 

Indeed, if p(x) > then (fTHll implies p(x) — u\p(x) \ — 7|„p(x) < f — \r\. If p(x) < 0, 
from (fE?|) it follows that r < p(x) +u\p(x)\ +l2nV( x )i hence — p(x) — u\p(x)\ — 7| n P'( a; ) — 

Next, a small computation proves that 

cond(p, x) < 7Tp^ u 72n~ 2 if and onl y if iLp( x ) < ^ [i 1 ~ u )b(^)l - tLp( x )] ■ 
Finally, from (JHJ) and dHJ) it follows 

\c-c\ <-i 2 2n p{x) < I [(l-u)|p(x)| -ll n p{x)] <i^\r\- 

From Lemma El we deduce that r is a faithful rounding of p(x). □ 

Numerical values of condition numbers for a faithful polynomial evaluation in IEEE- 
754 double precision are presented in Tabled for degrees varying from 10 to 500. 
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Table 1: A priori bounds on the condition number to ensure faithful rounding in IEEE- 
754 double precision for polynomials of degree 10 to 500 

n 10 100 200 300 400 500 

l^u7 2r r 2 1.13 • 10 13 1.13 ■ 10 11 2.82 • 10 1U 1.13 ■ 10 10 7.04 • 10 9 4.51 ■ 10 9 



4 Dynamic and validated error bounds for faithful 
rounding and accuracy 

The results presented in Section|3]are perfectly suited for theoretical purpose, for instance 
when we can a priori bound the condition number of the evaluation. However, neither 
the error bound in Theorem nor the criterion proposed in Theorem can be easily 
checked using only floating point arithmetic. Here we provide dynamic counterparts of 
Theorem 0] and Proposition [TJ that can be evaluated using floating point arithmetic in 
the "round to the nearest" rounding mode. 

Lemma 8. Consider a polynomial p of degree n with floating point coefficients, and x a 
floating point value. We use the notations of Algorithm® and we denote (p n + p a )(x) by 
c. Then 

c — c] < fl ; — := a. 18 

V l-2(n + l)u J y 1 

Proof. Let us denote Horner (\pn\ © \p<r\, \%\) by b. Since c = (p n + p a )(x) and c = 
Horner (p n © p a , x) where p n and p a are two polynomials of degree n—1, Lemma [T] yields 

\C~C\< 72n-l(& + Va){x) < (1 + u) 2n " 1 72n-l 6. 

From (HJ) and (JHJ) it follows that 

\C- C\ < (1 +U) 2 "72n-1& < (l + Ur+^^fo). 

Finally we use relation (j3J) to obtain the error bound. □ 

Remark 2. Lemma El allows us to compute a validated error bound for the computed 
correcting term c. We apply this result twice to derive next Theorem El First with 
Lemma El it yields the expected dynamic condition for faithful rounding. Then from the 
EFT for the Horner algorithm (Theorem EJ) we know that p(x) = P + c. Since r = r © c, 
we deduce \f — p(x)\ = \(r® c) — (r+ c) + (c — c)\. Hence we have 

\f-p(x)\ < |(f © c) - (f + c)| + |(c- c)|. (19) 

The first term |(r © c) — (r + c)| in the previous inequality is basically the absolute 
rounding error that occurs when computing r = r © c. Using only the bound (JHJ) of the 
standard model of floating point arithmetic, it could be bounded by u| r|. But here we 
benefit again from error free transformations using algorithm TwoSum to compute the 
actual rounding error exactly, which leads to a sharper error bound. Next Relation (120)) 
improves the dynamic bound presented in [2j. 
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Theorem 9. Consider a polynomial p of degree n with floating point coefficients, and x 
a floating point value. Let f be the computed value, f = CompHorner (p, x) (Algorithm^ 
and let a be the error bound defined by Relation Mt^l . 

i) If a < ^\r\, then r is a faithful rounding of p(x) . 

ii) Let e be the floating point value such that r + e = r + c, i.e., [f, e] = TwoSum (r, c), 
where r and c are defined by Algorithm^ The absolute error of the computed result 
f = CompHorner (p, x) is bounded as follows, 

\r- p ( x )\<fi(2±hy.= p. (20) 

Proof. The first proposition follows directly from Lemma |3 

By hypothesis r = ^r+c — e, and from Theorem |3] we have p(x) = r + c, thus 

| f — p(x)\ = | c — c — e\ < | c — c\ + \e\ < a + \e\. 

From (JHJ) and (jHJ) it follows that 

\f-p(x)\ < (1 + u)fl(S+ |e|) < fl 

which proves the second proposition. □ 

From Theorem |H1 we deduce the following algorithm. It computes the compensated 
result r together with the validated error bound (3. Moreover, the boolean value isfaithfu I 
is set to true if and only if the result is proved to be faithfully rounded. 

Algorithm 7. Compensated Horner algorithm with check of the faithful rounding 

function [r, (3, isfaithfu I] = CompHornerlsFaithul (p, x) 
[r,Pn,Pcr] = EFTHorner (p,x) 
c = Horner (p w © p a , x) 
b = Horner (\p w \ © \p a \, \x\) 
[r, e] = TwoSum ( r, c) 

a = (7an-i ® 6) (1 © 2(n + 1) ® u) 
P= (a© |e|) (1 - 2©u) 
isfaithful = (a<j\r\) 

5 Experimental results 

We consider polynomials p with floating point coefficients and floating point entries x. For 
presented accuracy tests we use Matlab codes for CompHorner (Algorithm lUJ) and Com- 
pHornerlsFaithul (Algorithm EJ) . These Matlab programs are presented in Appendix [7| 
From these Matlab codes, we see that CompHorner requires 0(21n) flop and that Com- 
pHornerlsFaithul requires 0(26n) flop. 

For time performance tests previous algorithms are coded in C language and several 
test platforms are described in next Table El 



/ a + |e|\ 
\l-2u) 
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Evaluation of p n (x) = (1-x) n in expanded form 

1.4e-11 
1.2e-11 

1e-11 

8e-12 

6e-12 

4e-12 

2e-12 



0.9 0.95 1 1.05 1.1 

1e+30 
1e+25 
% 1e+20 

Q. 

I 1e+15 

o 

1e+10 
100000 

0.9 0.95 1 1.05 1.1 

argument x 

Figure 1: We report the evaluation of polynomials p n near the multiple root x = 1 with the 
compensated Horner algorithm (CompHornerlsFaithul) and for multiplicity n — 6,8, 10, 12. 
Each evaluation proved to be faithfully rounded thanks to the dynamic test is reported 
with a green cross. The faithful evaluations that are not detected to be so with the 
dynamic test are represented in blue. Finally, the evaluations that are not faithfully 
rounded are reported in red. The lower frame represents the condition number with 
respect to the argument x. 

5.1 Accuracy tests 

We start testing the efficiency of faithful rounding with compensated Horner algorithm 
and the dynamic control of faithfulness. Then we focus more on both the a priori and 
dynamic bounds with two other test sets. Three cases may occur when the dynamic test 
for faithful rounding in Algorithm is performed. 

1. The computed result is faithfully rounded and this is ensured by the dynamic test. 
Corresponding plots are green in next figures. 

2. The computed result is actually faithfully rounded but the dynamic test fails to 
ensure this property. Corresponding plots are blue. 

3. The computed result is not faithfully rounded and plotted in red in this case. 
Next figures should be observed in color. 

5.1.1 Faithful rounding with compensated Horner 

In the first experiment set, we evaluate the expanded form of polynomials p n {x) = (1— x) n , 
for degree n — 6,8, 10, 12, at 2048 equally spaced floating point entries being near the 














n=12 






n=10 






n=8 






n=6 



12 



Accuracy of polynomial evaluation with the compensated Horner scheme [n=50] 




100000 1e+10 1e+15 1e+20 1e+25 1e+30 1e+35 

condition number 



Figure 2: We report the relative accuracy of every polynomial evaluation (y axis) with 
respect to the condition number (x axis). Evaluation is performed with CompHornerls- 
Faithul (Algorithm[7|). The color code is the same as for Figure^ Leftmost vertical line is 
the a priori sufficient condition ()16)l while the right one marks the inverse of the working 
precision u. Broken line is the a priori accuracy bound (JEJ). 



multiple root x = 1. These evaluations are extremely ill-conditioned since 

cond(j9 n ,x) = 

These condition numbers are plotted in the lower frame of Figure Q while x varies 
around the root. These huge values have a sense since polynomials p are exact in 
IEEE-754 double precision. Results are reported on Figure ^ The well known relation 
between the lost of accuracy and the nearness and the multiplicity of the root, i.e., the 
increasing of the condition number, is clearly illustrated. These results also illustrate 
that the dynamic bound becomes more pessimistic as the condition number increases. 
In next figures the horizontal axis does not represent the x entry range anymore but the 
condition number which governs the whole behavior. 

For the next experiment set, we first designed a generator of arbitrary ill-conditioned 
polynomial evaluations. It relies on the condition number definition ©. Given a degree n, 
a floating point argument x and a targeted condition number C, it generates a polynomial 
p with floating point coefficients such that cond(p, x) has the same order of magnitude 
as C. The principle of the generator is the following. 

1. \n/2\ coefficients are randomly selected and generated such that p(x) = 

l Oil i^i* ~ c , 

2. the remaining coefficients are generated ensuring \p{x)\ ~ 1 thanks to high accuracy 
computation. 

Therefore we obtain polynomials p such that cond(j», x) = p(x)/\p(x)\ ~ C, for arbitrary 
values of C. 
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Accuracy of the absolute error bounds for CompHorner 

1e-25 
1e-26 
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™ 1e-31 
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1e-34 

0.994 0.996 0.998 1 1.002 1.004 1.006 

argument x 

Figure 3: The dynamic error bound ()20|) compared to the a priori bound (|T3"|) and to the 
actual forward error (p(x) = (1 — x) 5 for 400 entries on the x axis). 

In this test set we consider generated polynomials of degree 50 whose condition num- 
bers vary from about 10 2 to 10 35 . These huge condition numbers again have a sense here 
since the coefficients and the argument of every polynomial are floating point numbers. 
The results of the tests performed with CompHornerlsFaithul (Algorithm |7|) are reported 
on Figure El As expected every polynomial with a condition number smaller than the a 
priori bound (|T^|) is faithfully evaluated with Algorithm [7| — green plots at the left of the 
leftmost vertical line. 

On Figure El we also see that evaluations with faithful rounding appear for condition 
numbers larger than the a priori bound (fTo^) — green and blue plots at the right of 
the leftmost vertical line. As expected a large part of these cases are detected by the 
dynamic test introduced in Theorem^ — the green ones. Next experiment set comes back 
to this point. We also notice that the compensated Horner algorithm produces accurate 
evaluations for condition numbers up to about 1/u — green and blue plots. 

5.1.2 Significance of the dynamic error bound 

We illustrate the significance of the dynamic error bound ()20|). compared to the a priori 
error bound (|T3"|) and to the actual forward error. We evaluate the expanded form of 
p(x) = (1 —x) 5 for 400 points near x — 1. For each value of the argument x, we compute 
CompHorner (p, x) (Algorithm |HJ), the associated dynamic error bound ()20j) and the actual 
forward error. The results are reported on Figure El 

As already noticed, the closer the argument is to the root 1 (i.e., , the more the 
condition number increases), the more pessimistic becomes the a priori error bound. 
Nevertheless our dynamic error bound is more significant than the a priori error bound 
as it takes into account the rounding errors that occur during the computation. 
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Table 2: Measured time performances for CompHorner, CompHornerlsFaithul and 
DDHorner. GCC denotes the GNU Compiler Collection and ICC denotes the Intel C/C++ 
Compiler. 





CompHorner 


CompHornerlsFaith 


DDHorner 


Horner 


Horner 


Horner 


Pentium 4, 3.00 GHz GCC 3.3.5 

ICC 9.1 


3.77 
3.06 


5.52 
5.31 


10.00 

8.88 


Athlon 64, 2.00 GHz GCC 4.0.1 


3.89 


4.43 


10.48 


Itanium 2, 1.4 GHz GCC 3.4.6 

ICC 9.1 


3.64 
1.87 


4.59 
2.30 


5.50 
8.78 




~ 2-4 


~4-6 


~ 5-10 



5.2 Time performances 

All experiments are performed using IEEE-754 double precision. Since the double- 
doubles j31 Ej are usually considered as the most efficient portable library to double the 
IEEE-754 double precision, we consider it as a reference in the following comparisons. 
For our purpose, it suffices to know that a double-double number a is the pair (a^, a{) of 
IEEE-754 floating point numbers with a = a^ + ai and |aj| < u|a/j|. This property implies 
a renormalisation step after every arithmetic operation with double-double values. We 
denote by DDHorner our implementation of the Horner algorithm with the double-double 
format, derived from the implementation proposed in 

We implement the three algorithms CompHorner, CompHornerlsFaith and DDHorner 
in a C code to measure their overhead compared to the Horner algorithm. We program 
these tests straightforwardly with no other optimization than the ones performed by the 
compiler. All timings are done with the cache warmed to minimize the memory traffic 
over- cost. 

We test the running times of these algorithms for different architectures with different 
compilers as described in Table 121 Our measures are performed with polynomials whose 
degree vary from 5 to 200 by step of 5. For each algorithm, we measure the ratio of its 
computing time over the computing time of the classic Horner algorithm; we display the 
average time ratio over all test cases in Table El 

The results presented in Table 121 show that the slowdown factor introduced by Com- 
pHorner compared to the classic Horner roughly varies between 2 and 4. The same slow- 
down factor varies between 4 and 6 for CompHornerlsFaithul and between 5 and 10 for 
DDHorner. We can see that CompHornerlsFaithul runs a most 2 times slower than Com- 
pHorner: the over-cost due to the dynamic test for faithful rounding is therefore quite 
reasonable. Anyway CompHorner and CompHornerlsFaithul run both significantly faster 
than DDHorner. 

Remark 3. We provide time ratios for IA'64 architecture (Itanium 2). Tested algorithms 
take benefit from IA'64 instructions, e.g., fma, but are not described in this paper. 
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6 Conclusion 



Compensated Horner algorithm yields more accurate polynomial evaluation than the clas- 
sic Horner iteration. Its accuracy behavior is similar to an Horner iteration performed in 
a doubled working precision. Hence compensated Horner may perform a faithful poly- 
nomial evaluation with IEEE- 754 floating point arithmetic in the "round to the nearest" 
rounding mode. An a priori sufficient condition with respect on the condition number 
that ensures such faithfulness has been defined thanks to the error free transformations. 

These error free transformations also allow us to derive a dynamic sufficient condition 
that is more significant to check for faithful rounding with compensated Horner algorithm. 

It is interesting to remark here that the significance of this dynamic bound can be 
improved easily — how to transform blue plots in green ones? Whereas bounding the 
error in the computation of the (polynomial) correcting term in Relation (|18j) . a good 
approximate of the actual error could be computed (applying again CompHorner to the 
correcting term). Of course such extra computation will introduce more running time 
overhead not necessary useful — green plots are here! So it suffices to run such extra (but 
costly) checking only if the previous dynamic one fails (a similar strategy as in dynamic 
filters for geometric algorithms). 

Compared to the classic Horner algorithm, experimental results exhibit reasonable 
over-costs for accurate polynomial evaluation (between 2 and 4) and even for this com- 
putation with a dynamic checking for faithfulness (between 4 and 6). Let us finally 
remark than such computation that provides as accuracy as if the working precision is 
doubled and a faithfulness checking is no more costly in term of running time than the 
"double-double" counterpart without any check. 

Future work will be to consider subnormals results and also an adaptative algorithm 
that ensure faithful rounding for polynomials with an arbitrary condition number. 
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7 Appendix 

Accuracy tests use next Matlab codes for algorithms Algorithm |H1 (CompHorner) and 
Algorithm [7| (CompHornerlsFaithul). Following Matlab convention, p is represented as a 
vector p such that p(x) = J2i=oV( n ~ * + We also recall that Matlab eps denotes 
the machine epsilon, which is the spacing between 1 and the next larger floating point 
number, hence u = eps/2. 
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Algorithm 8. Code for Algorithm |H1 

function r = CompHorner(p, x) 
n = length(p)-l; % degree of p 
[xh, xl] = Split(x); 
r = p(l) ; c = 0.0; 
for i=2:n+l 
%[r, pi] = TwoProd(r, x) 
p = r*x; 

[rh, rl] = Split(r); 

pi = rl*xl-(((p-rh*xl)-rl*xh)-rh*xl); 
%[r, sigma] = TwoSum(r, p(i)) 
r = p+p(i); 
t = r-p; 

sigma = (p-(r-t))+(p(i)-t); 
% Computation of the correcting term 
c = c*x+(pi+sig); 
end 

% Final correction of the result 
r = r+c; 



Algorithm 9. Code for Algorithm 

function [r, beta, isfaith] = CompHornerlsFaithul(p, x) 
n = length(p)-l; % degree of p 
[xh, xl] = Split(x); 
absx = abs(x); 

r = p (l); c = 0.0; beta = 0.0; 
for i=2:n+l 
% [ r , pi] = TwoProd(r, x) 

p = r*x; 

% [rh, rl] = Split(r); 

pi = rl*xl-(((p-rh*xl)-rl*xh)-rh*xl); 

% [r, sigma] = TwoSum(r, p(i)) 

r = p+p(i); 

t = r-p; 

sigma = (p-(r-t))+(p(i)-t); 
% Computation of the correcting term 
c = c*x+(pi+sig); 
b = b*absx+(abs(pi)+ abs(sig)); 
end 

% Final correction of the result 

[r, e] = TwoSum(r,c); 

% Check for faithful rounding 

alpha = gam(2*n-l)*b / (l-(n+l)*eps); 

isfaith = alpha j 0.25*eps*abs(r); 

% Absolute error bound 

beta = (alpha + abs(e))/(l-2*u); 
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